!
!  munilist
!
!  Copyright © 2011 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.


module mlist

  implicit none

  integer, parameter, private :: dbl = selected_real_kind(15)
  real(dbl), parameter, private :: rad = 57.295779513082322865_dbl


contains

  subroutine munilist(filename,output,column,idx,xcoo,ycoo,vxcoo,vycoo,tol, &
       min0,per,long,lat,alpha,delta,jd0,relative,mags,full,positions,header, &
       quiet,headkeys,fitskeys)

!    use cfitsio  ! consider fitsio
    use fitsio
    use trajd
    use traco

    character(len=*), intent(in) :: filename, output, column
    character(len=*), dimension(:), intent(in) :: fitskeys, headkeys
    integer, dimension(:), intent(in) :: idx
    real(dbl),dimension(:),intent(in) :: xcoo, ycoo, vxcoo, vycoo
    real(dbl),intent(in) :: tol, min0,per,long, lat, alpha, delta, jd0
    integer :: quiet
    logical, intent(in) :: mags,full,relative,positions,header

    integer, parameter :: DIM = 2
    integer :: status,blocksize,naxis,naxes(DIM),pcount,gcount,bitpix,nrows,&
         colnum
    character(len=10) :: colname
    logical :: simple,extend,anyf
    integer :: frow = 1, felem = 1
    real(dbl) :: nullval = 0.0

    character(len=80) :: filter, line, buf
    character(len=80), dimension(size(headkeys)) :: headval
    integer :: i,j,n,y1,m1,d1,h1,min1
    real(dbl) :: x,y,year,month,day,hour,min,sec,exptime,jd,ls,hcor,helcor, &
         phase,azimut,elev,ha,airm,r,rmin,dt
    real(dbl),dimension(:), allocatable :: flux,dflux,sky,dsky,xrect,yrect
    integer, dimension(:), allocatable :: id
    real(dbl), dimension(:), allocatable :: xpix,ypix,csky,cdsky,cflux,cdflux
!    real :: e,x(nstars),dx(nstars),sky(nstars),dsky(nstars)
!    real :: is(nstars),dis(nstars)

    n = size(idx) + size(xcoo)

    ! input FITS file
    status = 0
    call ftopen(20,filename,0,blocksize,status)
    call ftghpr(20,DIM,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
    if( status /= 0 ) then
       if( quiet == 0 ) write(*,*) "Failed to open: ", trim(filename)
       stop 66
    endif

    ! read date and time
    call ftgkys(20,fitskeys(1),line,buf,status)
    call fts2tm(line,y1,m1,d1,h1,min1,sec,status)
    year = y1; month = m1; day = d1; hour = h1; min = min1

    if( h1 == 0 .and. min1 == 0 ) then
       call ftgkys(20,fitskeys(2),line,buf,status)
       call fts2tm(line,y1,m1,d1,h1,min1,sec,status)
       hour = h1; min = min1
    endif

    ! read date, time, exposure, filter
    call ftgkye(20,fitskeys(3),exptime,buf,status)
    if( status /= 0 ) then
       exptime = 0.0
       if( quiet == 0 ) write(*,*) 'Exposure time not found. Assumed =0.'
       status = 0
    endif

    ! read date, time, exposure, filter
    call ftgkys(20,fitskeys(4),filter,buf,status)
    if( status /= 0 ) then
       filter = '-'
       if( quiet == 0 )  write(*,*) 'Filter not found.'
       status = 0
    endif

    ! required header quantities
    do i = 1, size(headkeys)

       call ftgkys(20,headkeys(i),headval(i),buf,status)
       if( status /= 0 ) then
          headval(i) = ''
          if( quiet == 0 )  write(*,*) "`",headkeys(i),"' not found."
          status = 0
       end if

    end do

    ! Evaluate Julian date at half of exposure 
    jd = datjd(year,month,day + (hour + (min + sec/60.0)/60.0)/24.0)
    jd = jd + exptime/2.0/86400.0
    
    ! select data extension
    call ftmnhd(20,BINARY_TBL,MEXTNAME,0,status)
    if( status == BAD_HDU_NUM ) then
       if( quiet == 0 )  write(*,*) 'MUNIPACK extension not found.'
    end if
    call ftgnrw(20,nrows,status)
!    call ftthdu(20,nhdu,stat)
!    do i = 1, nhdu
!       call ftmahd(20,i,hdutype,stat)
!       if( 

    ! read required data
    if( status == 0 .and. nrows > 0 ) then

       allocate(id(n))
       id = -1
       
       allocate(xpix(nrows),ypix(nrows))

       ! coordinates
       call ftgcnn(20,.false.,'x',colname,colnum,status)
       call ftgcvd(20,colnum,frow,felem,nrows,nullval,xpix,anyf,status)
       call ftgcnn(20,.false.,'y',colname,colnum,status)
       call ftgcvd(20,colnum,frow,felem,nrows,nullval,ypix,anyf,status)

       if( status /= 0 ) write(*,*) 'status=',status

       ! look for coordinates
       do j = 1, size(xcoo)
          dt = jd - jd0
          x = xcoo(j) + vxcoo(j)*dt
          y = ycoo(j) + vycoo(j)*dt
          rmin = huge(rmin)
          do i = 1, nrows
             r = sqrt((xpix(i) - x)**2 + (ypix(i) - y)**2)
             if( r < tol .and. r < rmin) then
                rmin = r
                id(j) = i
             end if
          end do

          if( id(j) < 0 ) write(*,*) 'Object at ',x,y,' not found.',rmin

       end do


       do j = 1, size(idx)
          id(j+size(xcoo)) = idx(j)
       end do

       allocate(csky(nrows),cdsky(nrows),cflux(nrows),cdflux(nrows))

       call ftgcnn(20,.false.,'sky',colname,colnum,status)
       call ftgcvd(20,colnum,frow,felem,nrows,nullval,csky,anyf,status)
       call ftgcnn(20,.false.,'sky_err',colname,colnum,status)
       call ftgcvd(20,colnum,frow,felem,nrows,nullval,cdsky,anyf,status)
       call ftgcnn(20,.false.,column,colname,colnum,status)
       call ftgcvd(20,colnum,frow,felem,nrows,nullval,cflux,anyf,status)
       call ftgcnn(20,.false.,trim(column)//'_err',colname,colnum,status)
       call ftgcvd(20,colnum,frow,felem,nrows,nullval,cdflux,anyf,status)

!       write(*,*) id

       allocate(sky(n),dsky(n),flux(n),dflux(n),xrect(n),yrect(n))
       do i = 1, n
          j = id(i)
          if( j > 0 ) then
             sky(i) = csky(j)
             dsky(i) = cdsky(j)
             flux(i) = cflux(j)
             dflux(i) = cdflux(j)
             xrect(i) = xpix(j)
             yrect(i) = ypix(j)
             !write(*,*) sky(i),dsky(i),flux(i),dflux(i)
          end if
       end do
       deallocate(csky,cdsky,cflux,cdflux)
       deallocate(id)
       deallocate(xpix,ypix)

    end if

    call ftclos(20,status)
    if( status /= 0 )then
       write(*,*) 'FITSIO: following error(s) has been occured:'
       do
          Call ftgmsg(buf)
          if( buf == ' ' ) exit
          write(*,*) trim(buf)
       enddo
       return
    endif



    ! heliocentric correction
    if( alpha >= 0.0 )then
       !  approx. of length of the Sun 
       ls = 279.465 + 0.985647*(jd - datjd(year,1.0_dbl,1.0_dbl))
       hcor = helcor(alpha/rad,delta/rad,ls/rad)
    else
       hcor = 0.0
    endif

    if( header ) then

       write(*,707) int(year),int(month),int(day),int(hour),int(min),int(sec),&
            jd,(trim(headval(i)),i=1,size(headval))
       
707    format (1x,i4.4,tr1,5(i2.2,tr1),f13.5,2000(1x,a))


    ! output relative quantities
    else if( relative ) then
       
       do i = 2, n
          if( flux(1) > 0.0 .and. flux(i) > 0.0 ) then
             dflux(i) = sqrt(dflux(i)**2 + dflux(1)**2)/flux(i) ! check!!!!
             flux(i) = flux(i)/flux(1)
          else
             flux(i) = -1.0
             dflux(i) = -1.0
          end if
       end do

       if( mags ) then
          where( flux > 0 )
             flux = 25.0 - 2.5*log10(flux)
             dflux = 1.086*dflux/flux
          elsewhere
             flux = 99.999
             dflux = 9.999
          end where
       end if

       write(*,710) jd + hcor,(flux(i),dflux(i),i=2,n)
710    format(1x,f13.5,1x,1000f7.3)

    else 
       ! instrumental fluxes
       
       if( full ) then

          if( alpha >= 0.0 .and. abs(lat) <= 90.0 )then
             ha = hangle(lmst(jd,long)*15.0_dbl,alpha)
             call eq2hor(ha,delta,lat,azimut,elev)
             airm = airmass(90.0_dbl-elev)
             write(*,701) int(year),int(month),int(day),int(hour),int(min), &
                  int(sec),exptime,filter,jd,azimut,elev,airm, &
                  (sky(i),dsky(i),flux(i),dflux(i),i = 1,n)
          else
!             write(*,*) n,sky,flux,dsky,dflux
             write(*,701) int(year),int(month),int(day),int(hour),int(min), &
                  int(sec),exptime,filter,jd,(sky(i),dsky(i),flux(i),dflux(i),i = 1,n)
          endif
701       format (1x,i4.4,tr1,5(i2.2,tr1),f6.1,tr1,a1,tr2,f13.5,1x, 1p, 1000g12.6)

       else if( positions .and. n > 0 ) then

          write(*,704) int(year),int(month),int(day),int(hour),int(min),int(sec), &
               jd, (xrect(i),yrect(i),i = 1,n)
704       format (1x,i4.4,tr1,5(i2.2,tr1),tr2,f13.5,1x, 10000f10.3)
          

       else

          write(*,700) int(year),int(month),int(day),int(hour),int(min),int(sec), &
               exptime,filter,jd, (flux(i),dflux(i),i = 1,n)
700       format (1x,i4.4,tr1,5(i2.2,tr1),f6.1,tr1,a1,tr2,f13.5,1x, 10000(es15.6,es9.2))


       end if



    end if

    deallocate(sky,dsky,flux,dflux)


  end subroutine munilist

end module mlist
